program PS

        ! "use" is the inheritance mechanism for F90:
        use mConstants
        use mObject
        use mUniverse
        use mLPS
        use mPSFittingFunctions
        use mMisc

        implicit none

        real(D) :: smithp3d, linearp3d, pdp3d
        real(D) :: k(200), kfactor, l1(200), smith1(200), pd1(200)
        real(D) :: sig8,dsig8,kr,Wsig,dlnk,twopisq
        real(D) :: om_m,om_v,om_b,h,A_s,n,w,wa
        integer :: i

        write( 6, '(a)' ) '# k(h/Mpc) P_Smith P_PD'

        kfactor = exp(log(1.e10_D)/199._D)
        k(1) = 1.e-5_D

        do i = 2, 200
          k(i) = k(i-1) * kfactor
        enddo
        
        om_m = 0.3_D
        om_v = 0.7_D
        om_b = 0.044_D
        h = 0.7_D
        A_s = 0.8_D
        n = 0.96_D
        w = -1.0_D
        wa = 0.0_D
        call def_universe(om_m,om_v,om_b,h,A_s,n,w,wa)

        twopisq = 19.7392_D
        do i = 1, 200
          smith1(i) = smithp3d(k(i),1._D);
          l1(i) = linearp3d(k(i),1._D);
          l1(i) = l1(i) * k(i)**3 /twopisq
          pd1(i) = pdp3d(k(i),1._D);
        enddo

        write( 6,*) '     k        Lin       Smith        PD'

        sig8 = 0._D
        do i = 1, 200
          kr = k(i)*(8.0_D)/h
          Wsig = 3._D*(sin(kr)-kr*cos(kr))
          dlnk = log(kfactor)

          dsig8 = W*W*l1(i)*dlnk
          sig8 = sig8 + dsig8
          write( 6, fmt = '(5e11.3)' ) k(i),l1(i),smith1(i),pd1(i), pd1(i)/smith1(i)
        end do
        write(6,*) 'sigma_8 = ',sig8

end program PS
